from time import time
import sys
sys.path.insert(0,'../..')
from rmpoly import *

def almost_equal(p1,p2):
  b = True
  assert p1.rp == p2.rp
  if len(p1) != len(p2):
    return False
  RR = p1.rp.field
  eps = RR(2^-RR.prec())
  for expv1,v1 in p1.iteritems():
    if abs(p2[expv1] - v1) > eps:
      return False
  return b

print '\ntest_poly_QQ.sage'

def test_str():
  rp,x,y = rgens('x,y',6,QQ)
  p = (1 + x + y)**2
  assert p == rp(str(p))

def test_gens():
  for ordx in ['lex','grlex',(0,0),(0,1),(1,0),(1,1)]:
    rp,x,y = rgens('x,y',6,QQ,order=ordx)
    assert x == rp('x')
    assert x**4 == rp('x^4')

def test_zero():
  rp = RPoly('x,y',6,QQ)
  p = rp('0')
  assert str(p) == '0'
  assert p == rp(0)

def test_code_expv():
  rp,x,y = rgens('x,y',6,QQ,order='grlex')
  t = (1,2)
  expv = rp.code_expv(t)
  assert rp.ngens == 2
  assert rp.mask_vars == (1<<14) - 1
  assert expv == 1 + (2<<7) + (3<<14)
  a = rp.decode_expv(expv)
  assert t == tuple(a)

  p = y^2
  expv = p.keys()[0]
  assert rp.decode_expv(expv) == [0,2]
  p = x^3
  expv = p.keys()[0]
  assert expv == 3 + (3<<14)
  # TODO is this the wanted behavior?
  #print rp.decode_expv(expv) #=> [3]
  #assert rp.decode_expv(expv) == [3,0]
  rp1,x,y = rgens('x,y',6,QQ)
  t = (1,2)
  expv = rp1.code_expv(t)
  assert rp1.ngens == 2
  assert expv == 1 + (2<<7)
  a = rp1.decode_expv(expv)
  assert t == tuple(a)

  p = y^2
  expv = p.keys()[0]
  assert rp.decode_expv(expv) == [0,2]

def test_new():
  rp = RPoly('x,y',6,QQ,order='grlex')
  rp1 = RPoly('x,y',9,QQ,order='grlex')
  p1 = rp('x^40*y + y^38')
  p2 = p1*p1
  assert p2.leading_expv()&rp.maskof != 0
  rp2 = RPoly(list('xy'),9,QQ,order='lex')
  p1a = rp2.new(p1)
  p1b = rp2('x^40*y + y^38')
  assert p1a == p1b
  p1 = rp1.new(p1)
  p2 = p1*p1
  assert p2.leading_expv()&rp1.maskof == 0

def test_mon_eval():
  rp = RPoly('x,y',6,QQ)
  p = rp.mon_eval('3*x*y^4')
  if rp.order == 'grlex':
    assert rp.mon_eval('3*x*y^4') == (int(1 + (4<<7) + (5<<14)), int(3))
  elif rp.order == 'lex':
    assert rp.mon_eval('3*x*y^4') == (int(1 + (4<<7)), int(3))


def test_from_mon():
  for ordx in ['lex','grlex',(0,0),(1,0),(0,1),(1,1)]:
    rp = RPoly('x,y',6,QQ,order=ordx)
    p = rp.from_mon((1,4,QQ('7/3')))
    assert p == rp('7/3*y^4')
    p = rp.from_mon(('y',4,QQ('7/3')))
    assert p == rp('7/3*y^4')

def test_from_expv():
  rp,x,y = rgens('x,y',6,QQ)
  p = x^2 + 2*y^2
  expv = p.leading_expv()
  p1 = rp.from_expv(expv)
  assert p1 == y^2

def test_variables():
   for ordx in ['lex','grlex', tuple([1]*11)]:
     rp = RPoly(['x%d' % i for i in range(11)],5,QQ,order=ordx)
     p = rp('72/5*x0^8*x1^6*x4^3*x10^2 + 2*x0^2*x2^3*x4^3')
     assert p.variables() == [0,1,2,4,10]

def test_lcm_expv():
  for ordx in ['lex','grlex',(0,0),(0,1),(1,0),(1,1)]:
    rp,x,y = rgens('x,y',6,QQ,order=ordx)
    p1 = x^2*y^4
    p2 = x^3*y
    expv1 = p1.keys()[0]
    expv2 = p2.keys()[0]
    expv = rp.lcm_expv(expv1,expv2)
    p = rp(0)
    p[expv] = 1
    assert p == x^3*y^4

def test_eq():
  rp,x,y = rgens('x,y',6,QQ)
  p1 = 1 + x
  p2 = 1 + 2*x
  p2 = p2 -x
  assert p1 == p2
  p2 -= x
  assert p2 == rp(1)
  assert p2 == int(1)
  assert p2 == 1

def test_add():
  for ordx in ['lex','grlex',(0,0),(0,1),(1,0),(1,1)]:
    rp,x,y = rgens('x,y',6,QQ,order=ordx)
    zero = rp(0)
    one = rp(1)
    assert zero + zero == zero
    assert zero + one == one
    p1 = -2*x
    p2 = 2*x
    assert p1+p2 == zero
    assert p1-p2 == -4*x
    p2 = 3/7*y^2
    p3 = p1 + p2
    assert p3 == -2*x+3/7*y^2
    p4 = p3
    p3 = p3 + 1
    assert p3 == -2*x+3/7*y^2 + 1
    assert p4 == -2*x+3/7*y^2
    p3 = p3 - 1
    assert p3 == -2*x+3/7*y^2
    p3 = 1 + p3
    assert p3 == -2*x+3/7*y^2 + 1
    p3 = 1 - p3
    assert p3 == 2*x-3/7*y^2
    p4 = x^2
    p4 = 1 - p4
    assert p4 == 1-x^2

  rp1,x0,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10 = rgens(['x%d' % i for i in range(11)],6,QQ)
  p1 = 36/5*x0^6*x1^6*x10^2 + x2^3
  p2 = 2*x0^2*x4^3
  assert p1+p2 == 36/5*x0^6*x1^6*x10^2 + 2*x0^2*x4^3 + x2^3

  rp2,x = rgens('x',8,QQ)
  p1 = -4*x^2 -2*x
  p2 = +2*x^2 +2*x +2
  assert p1 + p2 == -2*x^2 + 2

def test_iadd():
  rp,x,y = rgens('x,y',6,QQ)
  p3 = -2*x+3/7*y^2
  # check that p3 is mutable
  p4 = p3
  p3 += 1
  assert p3 == -2*x+3/7*y^2 + 1
  assert p4 == -2*x+3/7*y^2 + 1
  p3 = p3 - 1
  assert p3 == -2*x+3/7*y^2
  p3 -= 1
  assert p3 == -2*x+3/7*y^2-1
  p3 -= 1
  assert p3 == -2*x+3/7*y^2-2


def test_mul():
   rp,x = rgens('x',6,QQ)
   p1 = x^2
   assert p1 * 0 == 0
   assert 0*p1 == 0
   assert p1 * 1 == p1
   assert 1 * p1 == p1
   p2 = 2*x^60
   assert p1 * p2 == 2*x^62
   p1 = 1+x+x^2+x^3
   p2 = 1-x
   assert p1*p2 == 1-x^4

   rp1 = RPoly(['x%d' % i for i in range(11)],5,QQ)
   p1 = rp1('36/5*x0^6*x1^6*x10^2 + x2^3')
   p2 = rp1('2*x0^2*x4^3')
   p3 = p1*p2
   assert p3 == rp1('72/5*x0^8*x1^6*x4^3*x10^2 + 2*x0^2*x2^3*x4^3')


def test_pow_high():
  rp = RPoly(['x%d' % i for i in range(3)],6, QQ)
  x = rp.gens()
  p = 1 + sum(x)
  p1 = p**7
  p2 = p**14
  p3 = p.miller_pow(21)
  assert p3 == p1 * p2
  N = 10
  rp = RPoly(['x%d' % i for i in range(N)],6, QQ)
  x = rp.gens()
  p = 1 + sum(x)
  pw = 3
  p1 = p**pw
  p2 = p.miller_pow(pw)
  assert p1 == p2
  mc = multinomial_coefficients(len(p),pw)
  p3 = p.miller_pow(pw,mc)
  assert p2 == p3

def test_div():
  rp,x,y = rgens('x,y',6, QQ)
  p1 = x^2 + y^2 - 1
  p2 = p1/7
  assert p2 == 1/7*x^2 + 1/7*y^2 - 1/7

def test_mul_iadd():
  rp,x,y = rgens('x,y',6, QQ)
  p = x^2 + y^2 - 1
  p1 = x + y
  p2 = x - y
  p.mul_iadd(p1,p2)
  assert p == 2*x^2 - 1

def test_mul_mon():
  for ordx in ['lex','grlex',(0,0),(1,0),(0,1),(1,1)]:
    rp,x,y = rgens('x,y',6, QQ,order=ordx)
    p1 = x^2 + y^2 - 1
    p2 = p1.mul_mon((0,4,2/7))
    assert p2 == p1*2/7*x^4
    expv1 = (x^4).keys()[0]
    p2 = p1.mul_mon((expv1,2/7),1)
    assert p2 == p1*2/7*x^4

def test_iadd_mon():
  for ordx in ['lex','grlex',(0,0),(1,0),(0,1),(1,1)]:
    rp,x,y = rgens('x,y',6, QQ,order=ordx)
    p1 = x^2 + y^2
    p1.iadd_mon(((int(1),int(2)),2))
    p1a = x^2 + y^2 + 2*x*y^2
    assert p1 == p1a
    p1 = x^2 + y^2

    if rp.order == 'grlex':
      expv = int(1 + (2<<7) + (3<<14))
    elif rp.order == 'lex':
      expv = int(1 + (2<<7))
    if rp.order == 'grlex' or rp.order == 'lex': 
      p1.iadd_mon((expv,int(2)),int(1))
      assert p1 == p1a

def test_high_exp():
  rp,y,x = rgens('y,x',110,Integer)
  n = 2^1000
  p = x^n + y^n
  p1 = p**4
  assert p1.coeff(x^(2*n)*y^(2*n)) == 6

def test_total_degree():
  gens=['x%d' % i for i in range(11)]
  for ordx in ['lex','grlex',tuple([0]*9 + [1,1])]:
    rp = RPoly(gens,6,QQ,order=ordx)
    p = rp('36/5*x0^6*x1^6*x10^2 + 2*x0^2*x4^3 + x2^3')
    assert p.total_degree() == 14

def test_derivative():
  gens=['x%d' % i for i in range(11)]
  for ordx in ['lex','grlex',tuple([0]*11), tuple([1]*9 + [0,0])]:
    rp = RPoly(gens,6,QQ,order=ordx)
    p = rp('288/5*x0^8*x1^6*x4^3*x10^2 + 8*x0^2*x2^3*x4^3 +2*x0^2-2*x1^2')
    #print 'DB40', type(dict(p).items()[0][0])
    p1 = p.derivative(0)
    assert p1 == rp('2304/5*x0^7*x1^6*x4^3*x10^2 + 16*x0*x2^3*x4^3 + 4*x0')
    p1 = p.derivative('x0')
    assert p1 == rp('2304/5*x0^7*x1^6*x4^3*x10^2 + 16*x0*x2^3*x4^3 + 4*x0')

def test_integrate():
  for ordx in ['lex','grlex',(0,0),(0,1),(1,0),(1,1)]:
    rp,x,y = rgens('x,y',6,QQ,order=ordx)
    p = x^2*y^5 + x + y
    p1 = p.integrate('x')
    assert p1 == 1/3*x^3*y^5 +1/2*x^2 + x*y
    p1 = p.integrate('y')
    assert p1 == 1/6*x^2*y^6 + x*y + 1/2*y^2


def test_copy():
  for ordx in ['lex','grlex',(0,0),(0,1),(1,0),(1,1)]:
    rp = RPoly(['x','y'],6, QQ,order=ordx)
    p1 = rp('x^2 + y^2 - 1')
    p2 = p1.copy()
    p1 += int(1)
    assert p1 == rp('x^2 + y^2')
    assert p2 == rp('x^2 + y^2 - 1')

def test_coefficient():
  gens=['x%d' % i for i in range(11)]
  orders = ['lex','grlex',tuple([0]*11), tuple([0,0] + [1]*9)]
  for ordx in orders:
    rp = RPoly(gens,6,QQ,order=ordx)
    p1 = rp('x0^6 + 6*x0^5*x1 + 15*x0^4*x1^2 + 20*x0^3*x1^3 + 15*x0^2*x1^4 + 6*x0*x1^5 + x1^6 + x0^4*x1^2*x2')
    m = rp('x0^2')
    p2 = p1.coefficient(m)
    assert p2 == rp('15*x1^4')
    m = rp('x0^2*x1^2')
    p2 = p1.coefficient(m)
    assert p2 == rp('0')
    m = rp('x0^2*x1^4')
    p2 = p1.coefficient(m)
    assert p2 == rp('15')
    m = rp('x0^4*x1^2')
    p2 = p1.coefficient(m)
    assert p2 == rp('x2+15')
    assert p1.coefficient_t((1,1)) == p1.coefficient(rp('x1'))


def test_coeff():
  gens=['x%d' % i for i in range(11)]
  orders = ['lex','grlex',tuple([0]*11), tuple([0,0] + [1]*9)]
  for ordx in orders:
    rp = RPoly(gens,6,QQ,order=ordx)
    p1 = rp('x0^6 + 6*x0^5*x1 + 15*x0^4*x1^2 + 20*x0^3*x1^3 + 15*x0^2*x1^4 + 6*x0*x1^5 + x1^6')
    assert p1.coeff(rp('x0^2')) == int(0)
    assert p1.coeff(rp('x0^2*x1^4')) == int(15)


def test_iadd_m_mul_q():
  for ordx in ['lex','grlex',(0,0),(0,1),(1,0),(1,1)]:
    rp,x,y = rgens('x,y',6, QQ,order=ordx)
    p1 = x^2 + y^2
    p2 = x^3 + y^3
    m = x*y
    expva = m.keys()[0]
    c = m[expva]
    p1.iadd_m_mul_q(p2,(expva,c))
    assert p1 == x^2 + y^2 + x^4*y + x*y^4
    m = rp('x')
    expva = m.keys()[0]
    c = m[expva]
    p1.iadd_m_mul_q(p2,(expva,c))
    assert p1 == x^4*y + x*y^4 + x^4 + x*y^3 + x^2 + y^2
    m = -x*y
    expva = m.keys()[0]
    c = m[expva]
    p1.iadd_m_mul_q(p2,(expva,c))
    assert p1 == x^4 + x*y^3 + x^2 + y^2

def test_mul_trunc():
  for ordx in ['lex','grlex',(0,0,0),(0,1,0),(1,0,0),(1,1,0)]:
    rp,t,y,x = rgens('t,y,x',6,QQ,order=ordx)
    p1 = 1 + t^6*x + t^2*y
    p2 = p1.mul_trunc(p1,'t',5)
    assert p2 == y^2*t^4 + 2*y*t^2 + 1
    p2 = p1*p1
    p2 = p2.trunc('t',5)
    assert p2 == y^2*t^4 + 2*y*t^2 + 1
    p = 1 + t*x + t*y
    for i in range(2):
      p = p.mul_trunc(p,'t',3)
    assert p == 6*x^2*t^2 + 12*x*y*t^2 + 6*y^2*t^2 + 4*x*t + 4*y*t + 1
  
  for ordx in ['grlex',(1,1,1)]:
    rp,t,y,x = rgens('t,y,x',6,QQ,order=ordx)
    p = 6*x^2*t^2 + 12*x*y*t^2 + 6*y^2*t^2 + 4*x*t + 4*y*t + 1
    for i in range(2):
      p = p.mul_trunc(p,'_t',6)
    assert p == 120*x^2*t^2 +240*x*y*t^2 +120*y^2*t^2 +16*x*t +16*y*t +1
    p = 1 + x + y^2 + t^3
    p1 = rp(1)
    for i in range(10):
      p1 = p1.mul_trunc(p,'_t',6)
    assert p1 == 252*x^5 + 840*x^3*y^2 + 360*x*y^4 + 360*x^2*t^3 + 90*y^2*t^3 + 210*x^4 + 360*x^2*y^2 + 45*y^4 + 90*x*t^3 + 120*x^3 + 90*x*y^2 + 10*t^3 + 45*x^2 + 10*y^2 + 10*x + 1

  rp,a,b,y,x = rgens('a,b,y,x',6,QQ,order=(1,1,0,0))
  p = 1 + a*x*y + b*y
  p1 = p
  for i in range(2):
    p1 = p1.mul_trunc(p1,'_t',3)
  assert p1 == 6*x^2*y^2*a^2 +12*x*y^2*b*a +6*y^2*b^2 +4*x*y*a +4*y*b +1

  rp,a,b,y,x = rgens('a,b,y,x',8,QQ,order=(1,2,0,0))
  p = 1 + a*x*y + b*y
  p1 = p
  for i in range(2):
    p1 = p1.mul_trunc(p1,'_t',3)
  assert p1 == 6*x^2*y^2*a^2 +4*y*b +4*x*y*a +1
  assert p.pow_trunc(4,'_t',3) == p1
  assert p.pow_trunc(100,'_t',3) == \
      4950*x^2*y^2*a^2 +100*y*b +100*x*y*a +1

  rp,y,x = rgens('y,x',6,QQ,order=ordx)
  p1 = 1 + x + y
  p2 = p1.pow_trunc(3,'y',1)
  assert p2 == x^3 +3*x^2 +3*x +1
  p2 = p1.pow_trunc(4,'y',2)
  assert p2 == (p1^4).trunc('y',2)
  rp,y,z,x = rgens('y,z,x', 10,QQ)
  p3 = 1+x+x^2+x^10+x^20
  p = p3.pow_trunc(1,'x',4)
  assert p == 1+x+x^2

  for ordx in ['lex','grlex',(0,0,0),(0,1,0),(1,0,0),(1,1,0)]:
    rp,t,y,x = rgens('t,y,x',10,QQ,order=ordx)
    p1 = 1 + t^6*x + t^2*y
    p2 = p1.mul_trunc(p1,['t'],[5])
    assert p2 == y^2*t^4 + 2*y*t^2 + 1
    p2 = p1.mul_trunc(p1,['t','x'],[5,3])
    assert p2 == y^2*t^4 + 2*y*t^2 + 1
    p2 = p1.mul_trunc(p1,['t','x','y'],[5,3,2])
    assert p2 == 2*y*t^2 + 1
    p1 = p1**2
    p2 = p1*p1
    for hv in [[3,3,3],[3,4,5],[5,4,3],[4,5,6]]:
      p2b = p1.mul_trunc(p1,['t','x','y'],hv)
      p2a = p1*p1
      p2a = p2a.trunc('t',hv[0])
      p2a = p2a.trunc('y',hv[1])
      p2a = p2a.trunc('x',hv[2])
      assert p2b == p2a

  rp,x,y = rgens('x,y',6,QQ)
  p1 = (1 + x + x*y + 2*x*y^3)^3
  p2 = (1 + x + x*y + 3*x*y^2)^3
  p3 = p1.mul_trunc(p2,['x','y'],[2,7])
  p4 = p1.mul_trunc(p2,['y','x'],[7,2])
  assert p3 == p4


def test_trunc():
  rp = RPoly(['y','x'],6,QQ)
  p1 = rp('1 + x + y')**6
  p2 = p1.trunc('x',4).trunc('y',4)
  p3 = p1.trunc(['x','y'],[4,4])
  assert p2 == p3

def test_Subs():
  gens0=['x%d' % i for i in range(11)]
  orders = ['lex','grlex',tuple([0]*11), tuple([0,0] + [1]*9)]
  #orders = ['lex','grlex']
  for ordx in orders:
    rp = RPoly(gens0,6,QQ,order=ordx)
    gens = rp.gens()
    p1 = rp('x0^2*x1 + 2*x0^2 + 3*x1+ 7*x3*x4*x5 + x6*x7*x8')
    rp1 = RPoly(['c'],6,QQ)
    c = rp1.gens()[0]
    rules = {'x0':c,'x1':c+1,'x2':c^2+1,'x3':c+2,'x4':c^4,
    'x5':c+1,'x6':c*(c-3),'x7':c*(c-7),'x8':c+7,
    'x9':c+9, 'x10':c+10}
    sb = Subs(rp, rp1,rules)
    p = sb.subs(p1)
    assert p == 7*c^6 + 22*c^5 + 11*c^4 - 48*c^3 + 150*c^2 + 3*c + 3
  
    p = sb.subs_trunc(p1,'c',3)
    assert p == 150*c^2 +3*c +3

    rp,z,y,x = rgens('z,y,x',6,QQ)
    p1 = x^2 + 2*y^2 + 3*z^2
    rules = {'z': x+y}
    sb = Subs(rp, rp, rules)
    p2 = sb.subs(p1)
    p3 = x^2 + 2*y^2 + 3*(x+y)^2
    assert p2 == p3

  rp,y,x = rgens('y,x',6,QQ)
  rules = {'x':(1+y)^5}
  sb = Subs(rp,rp,rules)
  p = (1+x)^4
  p1 = sb.subs(p)
  assert p1.coefficient(y^20) == rp(1)
  p2 = sb.subs_trunc(p,'y',4)

  rules = {'x':(1+y)^5}
  sb = Subs(rp,rp,rules)
  p = (1+x)^4
  p2 = sb.subs_trunc(p,'y',4)
  p1 = sb.subs(p)
  assert p1.coefficient(y^20) == rp(1)

  sb = Subs(rp,rp,{'x':1+y^2,'y':1+x})
  p = 1 + x^2 + 2*y^2
  p1 = sb.subs(p)
  assert p1 == 1 + (1+y^2)^2 + 2*(1+x)^2

  # conversion from polynomial on QQ and polynomial on RealField(prec)
  RR = RealField(100)
  rp,y,x = rgens('y,x',6,QQ)
  rp1,c = rgens('c',6,RR)
  rules = {'x':RR(1/3)*c, 'y':RR(1/7)}
  sb = Subs(rp, rp1, rules)
  p = (x+y)^6
  p1 = sb.subs(p)
  p2 = (RR(1/3)*c + RR(1/7))**6
  assert p1.coeff(c) == p2.coeff(c)
  assert almost_equal(p1,p2)


def test_subs():
  rp,z,y,x = rgens('z,y,x',6,QQ)
  p1 = x^2 + 2*y^2 + 3*z^2
  p2 = p1.subs(z=x+y)
  p3 = x^2 + 2*y^2 + 3*(x+y)^2
  assert p2 == p3

  rp,y,x = rgens('y,x',6,QQ)
  p = (1+x)^4
  p1 = p.subs(x=(1+y)^5)
  assert p1.coefficient(y^20) == rp(1)

  p = 1 + x**2 + 2*y**2
  p1 = p.subs(x=1+y**2,y=1+x)
  assert p1 == 1 + (1+y**2)**2 + 2*(1+x)**2


def test_leading_expv():
  rp,x,y = rgens('x,y',6, QQ)
  p = 2*x^2 + 3*y^2 + x + y + 1
  expv = p.leading_expv()
  assert rp.decode_expv(p.leading_expv()) == [0,2]
  p1 = rp(0)
  p1[expv] = 1
  assert p1 == y^2
  p = 2*x^3 + 3*y^2 + x + y + 1
  expv = p.leading_expv()
  if rp.order == 'grlex':
    assert expv == int(3 + (3<<14))
  elif rp.order == 'lex':
    assert expv == int(2<<7)
  rp1,x,y,z,w = rgens('x,y,z,w',6,QQ)
  p = x^2*y*z^2 + 2*x^2*y*z^2 + 3*x^2*y^2*z
  expv = p.leading_expv()
  assert p[expv] == 3
  p = x^3*y*z^2 + 2*x^2*y*z^2 + 3*x^2*y^2*z
  expv = p.leading_expv()
  assert p[expv] == 1
  for ordx in ['lex',(0,0,0),(0,0,1)]:
    rp1,z,y,x = rgens('z,y,x',6,Integer,order=ordx)
    p1 = x^3*y^2*z
    p2 = x*y^2*z^4
    assert p1.leading_expv() > p2.leading_expv()
    assert x.leading_expv() > (y^5).leading_expv() 

  for ordx in ['grlex',(1,1,1),(1,1,0),(1,0,1)]:
    rp2,z,y,x = rgens('z,y,x',6,Integer,order=ordx)
    p1 = x^3*y^2*z
    p2 = x*y^2*z^4
    assert p2.leading_expv() > p1.leading_expv()

 
def test_leading_term():
  for ordx in ['lex','grlex',(0,0),(0,1),(1,1)]:
    rp,x,y = rgens('x,y',6, QQ)
    p = 2*x^2 + 3*y^2 + x + y
    p1 = p.leading_term()
    assert p1 == 3*y^2

  rp1,z,y,x = rgens('z,y,x',6, QQ)
  p1 = 2*x^2*y^8 - 3*x^5*y*z^4 + x*y*z^3 - x*y^4
  assert rp1.decode_expv(p1.leading_expv()) == [4,1,5]

  for ordx in ['lex','grlex',(0,0,0),(0,0,1),(1,0,1),(1,1,1)]:
    rp2,z,y,x = rgens('z,y,x',6,QQ,order=ordx)
    p2 = 2*x^2*y^8 - 3*x^5*y*z^4 + x*y*z^3 - x*y^4
    assert p2.leading_term() == -3*x^5*y*z^4

def test_division():
  for ordx in ['lex',(0,0),(0,1),(1,2)]:
    rp,y,x = rgens('y,x',6, QQ,order=ordx)
    f = x^3
    f0 = x-y^2
    f1 = x-y
    qv, r = f.division((f0,f1))
    assert (qv[0],qv[1],r) == (x^2 + x*y^2 + y^4,rp(0),y^6)
  for ordx in ['grlex',(1,1),(2,3)]:
    rp,y,x = rgens('y,x',6, QQ,order=ordx)
    f = x^3
    f0 = x-y^2
    f1 = x-y
    qv, r = f.division((f0,f1))
    assert (qv[0],qv[1],r) == (-x-1,x^2+x*y+x+y+1,y)
  rp,y,x = rgens('y,x',6, QQ,order=(1,0))
  f = x^3
  f0 = x-y^2
  f1 = x-y
  qv, r = f.division((f0,f1))
  assert (qv[0],qv[1],r) == (rp(0),rp(0),x^3)

  for ordx in ['lex',(0,0)]:
    rp1,x,y = rgens('x,y',6, QQ, order=ordx)
    f = x^3*y^2 + x*y^4
    g = x^2 + y
    qv, r = f.division((g,))
    assert r == x^9 +x^7
    f = x^10*y^2 + x^20*y^4
    qv, r = f.division((g,))
    assert r == x^28 + x^14
    f = x^20*y^2 + 2*x^25*y^4
    qv, r = f.division((g,))
    assert r == 2*x^33 + x^24
    f = x^61*y^2 + x^60*y^4
    bof = 0
    try:
      qv, r = f.division((g,))
    except RPolyOverflowError:
      bof = 1
    assert bof == 1
    rp2,x,y = rgens('x,y',7, QQ)
    f = x^61*y^2 + x^60*y^4
    g = x^2 + y
    qv, r = f.division((g,))
    assert r == x^68 + x^65

def test_reversion():
  rp,y,x = rgens('y,x',6,QQ)
  p = x + x**2
  n = 10
  r = p.series_reversion('x',n, 'y')
  assert r == 1430*y^9 -429*y^8 +132*y^7 -42*y^6 +14*y^5 -5*y^4 +2*y^3 -y^2 +y
  p = x + 2*x^2 + x^3
  r = p.series_reversion('x',n, 'y')
  assert r == y - 2*y^2 + 7*y^3 - 30*y^4 + 143*y^5 - 728*y^6 + 3876*y^7 - 21318*y^8 + 120175*y^9
  p = 2*x + 2*x^2 + x^3
  r = p.series_reversion('x',n,'y')
  assert r == 1/2*y - 1/4*y^2 + 3/16*y^3 - 5/32*y^4 + 17/128*y^5 - 7/64*y^6 + 21/256*y^7 - 99/2048*y^8 + 55/8192*y^9

  rp,z,y,x = rgens('z,y,x',6,QQ)
  p = x + x**2 + x**2*y + x**3*y
  n = 10
  r = p.series_reversion('x',n, 'z')
  assert p.subs_trunc('z',n,x=r) == z
  assert r.subs(y=1) == z - 2*z^2 + 7*z^3 - 30*z^4 + 143*z^5 - 728*z^6 + 3876*z^7 - 21318*z^8 + 120175*z^9

def test_inversion():
  rp,x = rgens('x',6,QQ)
  p = 2 + x + 2*x^2
  n = 5
  p1 = p.series_inversion('x',n)
  #print 'DB42', p*p1, (p*p1).trunc('x',n)
  assert (p*p1).trunc('x',n) == rp(1)
  rp,y,x = rgens('y,x',6,QQ)
  n = 5
  p = 2 + x*(1+y) + x^2*(2-3*y)
  p1 = p.series_inversion('x',n)
  assert (p*p1).trunc('x',n) == rp(1)

  rp,y,x = rgens('y,x',6,QQ,order='grlex')
  p = 1 + x + y
  n = 3
  p1 = p.series_inversion('_t',n)
  assert p1 == 1 - (x+y) + (x+y)^2

  rp,x,y = rgens('x,y',6,QQ)
  p = 2 + x + y + 2*(x+y)^2
  nx = 5
  ny = 5
  p1 = p.series_inversion(['x','y'],[nx,ny])
  rp2,x2 = rgens('x2',6,QQ)
  p2 = 2 + x2 + 2*x2^2
  n = nx+ny
  p2a = p2.series_inversion('x2',n)
  sb = Subs(rp2,rp,{'x2':x+y})
  p2a = sb.subs_trunc(p2a,['x','y'],[nx,ny])
  assert p1 == p2a

def test_reversion2():
  rp,z,y,x = rgens('z,y,x',10,QQ)
  p = x*(1+y) + x^2*(1+y)
  nx,ny = 4,4
  p1 = p.series_reversion2('x',nx,'y',ny,'z')
  assert p.subs_trunc(['z','y'],[nx,ny],x=p1) == z
  nx = 4
  ny = 5
  p = x*(3+y+y^2)
  for i in range(2,nx):
    p += x^i*(1+y)^i
  p1 = p.series_reversion2('x',nx,'y',ny,'z')
  assert p.subs_trunc(['z','y'],[nx,ny],x=p1) == z

  p = x*(1+y) + x^2*(1+y)
  nx,ny = 4,4
  p1 = p.series_reversion2('x',nx,'y',ny,'z')
  assert p.subs_trunc(['z','y'],[nx,ny],x=p1) == z

  rp,w,z,y,x = rgens('w,z,y,x',10,QQ)
  p = x*(1+y+z) + x^2*(1+y+z^2)
  nx,ny,nz = 4,4,4
  p1 = p.series_reversion2('x',nx,['y','z'],[ny,nz],'w')
  assert p.subs_trunc(['w','y','z'],[nx,ny,nz],x=p1) == w

def test_quaternions():
  R = QuaternionAlgebra(-1,-1)
  _,i,j,k = R.basis()
  rp,x1,x2,x3,y1,y2,y3 = ncrgens('x1,x2,x3,y1,y2,y3',6,R)
  v1 = x1*i + x2*j + x3*k
  v2 = y1*i + y2*j + y3*k
  assert v1^2 == -(x1^2 + x2^2 + x3^2)
  assert v1*v2 == -x1*y1-x2*y2-x3*y3 + i*(x2*y3-x3*y2) + j*(x3*y1-x1*y3) + \
      k*(x1*y2-x2*y1)

  assert v1*i == -x1 - x2*k + x3*j
  assert i*v1 == -x1 + x2*k - x3*j

  rp1,x = ncrgens('x',8,R)
  rp2,y = rgens('y',8,rp1)
  p = (x + i*x^2 + j*x^3)*y + k*x + 1
  # not working because
  # sage: Q.<i,j,k> = QuaternionAlgebra(-1,-1)
  # sage: Q('-1 + 2*j')  
  #p2 = rp2('+( +(-1)*x^6 +(-1 + 2*j)*x^4 +(2*i)*x^3 +(1)*x^2)*y^2 +( +(2*j)*x^3 +(2*i)*x^2 +(2)*x)*y +( +(1))')
  #assert p2 == p^2
  p2 = (-x^6 + (-1 + 2*j)*x^4 + 2*i*x^3 + x^2)*y^2 + (2*j*x^3 + 2*i*x^2 + 2*x + 2*k*x^2)*y + 1 - x^2 + 2*k*x
  assert p^2 == p2
  assert p.mul_trunc(p,'y',2) == (2*j*x^3 + 2*i*x^2 + 2*x + 2*k*x^2)*y + 1 - x^2+ 2*k*x
  p = i*x*y + j
  assert p^2 == -x^2*y^2 -1
  assert k*p == j*x*y - i
  p1 = x*k
  assert p1*p == j*x^2*y -i*x
 

def test_series():
  R.<x> = QQ[[]]
  p = 1 + x + O(x^6)
  rp,y = mrgens('y',8,R)
  p1 = p + p^2*y
  p2 = p1.pow_trunc(4,'y',4)
  assert p2 == (1 + 4*x + 6*x^2 + 4*x^3 + x^4 + O(x^6)) + (4 + 20*x + 40*x^2 + 40*x^3 + 20*x^4 + 4*x^5 + O(x^6))*y + (6 + 36*x + 90*x^2 + 120*x^3 + 90*x^4 + 36*x^5 + O(x^6))*y^2 + (4 + 28*x + 84*x^2 + 140*x^3 + 140*x^4 + 84*x^5 + O(x^6))*y^3
  assert str(p2) == ' +(4 + 28*x + 84*x^2 + 140*x^3 + 140*x^4 + 84*x^5 + O(x^6))*y^3 +(6 + 36*x + 90*x^2 + 120*x^3 + 90*x^4 + 36*x^5 + O(x^6))*y^2 +(4 + 20*x + 40*x^2 + 40*x^3 + 20*x^4 + 4*x^5 + O(x^6))*y +(1 + 4*x + 6*x^2 + 4*x^3 + x^4 + O(x^6))'



def test_matrix():
  M = MatrixSpace(QQ,2,2)
  rp,x = ncrgens('x',6,M)
  m1 = M([1,2,3,4])
  m2 = M([5,6,7,8])
  p = m1*x + m2
  assert p*m1 == m1^2*x + m2*m1
  assert m1*p == m1^2*x + m1*m2
  p2 = p^2
  assert p2.coeff(x^2) == m1^2
  assert p2.coeff(x) == m1*m2 + m2*m1
  assert p2[0] == m2^2
  assert p*m1 == m1^2*x + m2*m1
  

  rp,y,x = ncrgens('y,x',6,M)
  p = 1 + m1*x + m2*y
  p1 = p.log(['y','x'], [3,3])
  assert p1.coeff(x) == m1
  assert p1.coeff(y) == m2
  assert p1.coeff(x^2) == -1/2 * m1^2
  assert p1.coeff(x*y) == -1/2 * (m1*m2 + m2*m1)
  assert p1.coeff(x^2*y) == 1/3 * (m1^2*m2 + m1*m2*m1 + m2*m1^2)
  assert p1.coeff(x*y^2) == 1/3 * (m2^2*m1 + m2*m1*m2 + m1*m2^2)
  assert p1.coeff(x^2*y^2) == -1/4 * (m1^2*m2^2 + m2^2*m1^2 + (m1*m2 + m2*m1)^2)

  rp,x = ncrgens('x',6,M)
  p = 1 + m1*x + m2*x^2
  p2 = p.pow_trunc(2,'x',10)
  assert p2 == (1 + m1*x)^2 + 2*m2*x^2 + (m1*m2+m2*m1)*x^3 + m2^2*x^4
  p3 = p2.nth_root(2,'x',10)
  assert p3 == p
  p = 1 + m1*x + m2*x^2
  p2 = p.pow_trunc(3,'x',10)
  p3 = p2.nth_root(3,'x',10)
  assert p3 == p
  p1a = (x*m1).exp('x',4)
  # Baker Campbell Hausdorff formula
  h = 5
  p1 = (x*m1).exp('x',h).mul_trunc((x*m2).exp('x',h),'x',h)
  c2 = m1*m2-m2*m1
  c3 = m1*c2-c2*m1 - m2*c2+c2*m2
  c4 = m1^2*m2^2 - 2*m1*m2*m1*m2 - m2^2*m1^2 + 2*m2*m1*m2*m1
  p2 = x*(m1+m2) + x^2*c2/2 + x^3*c3/12 + x^4*c4/24
  p3 = p2.exp('x',h)
  assert p1 == p3
  p5 = p1.log('x',h)
  assert p5 == p2

  # series inversion
  p = m1 + x*m2
  p1 = p.series_inversion('x',3)
  n1 = m1^-1
  p2 = n1 - x*n1*m2*n1 + x^2*n1*m2*n1*m2*n1
  assert  p1 == p2


  p = m1*x + m2*x^2
  p1 = p.atan('x',h)
  p2 = p1.tan('x',h)
  assert p == p2
  p3 = p1.cos('x',h)
  assert p1.sin('x',h) == p2.mul_trunc(p3,'x',h)


  p1 = p.atanh('x',h)
  p2 = p1.tanh('x',h)
  assert p == p2
  p3 = p1.cosh('x',h)
  assert p1.sinh('x',h) == p2.mul_trunc(p3,'x',h)
  p1 = p.atanh('x',h)
  p2 = p1.tanh('x',h)
  assert p == p2

  M = MatrixSpace(QQ,2,2)
  rp,x = ncrgens('x',6,M)
  m1 = M([1,2,3,4])
  p = 1 + m1*x
  p1 = p.log('x',3)
  assert p1.coeff(x) == m1
  assert p1.coeff(x^2) == -1/2*m1^2

def test_free_algebra():
  F.<A,B> = FreeAlgebra(QQ,2)
  rp,x,y = ncrgens('x,y',10,F, order='grlex',base_field=QQ)
  # test with the Baker-Campbell-Hausdorff formula
  # http://en.wikipedia.org/wiki/Baker%E2%80%93Campbell%E2%80%93Hausdorff_formula
  h = 5
  pa = (x*A).exp('_t',h)
  pb = (y*B).exp('_t',h)
  p = pa.mul_trunc(pb,'_t',h)
  p1 = p.log('_t',h)
  cxy = 1/2*(A*B - B*A)
  cxy2 = 1/12*(A*B^2 - 2*B*A*B + B^2*A)
  cx2y = 1/12*(A^2*B - 2*A*B*A + B*A^2)
  cx2y2 = 1/24*(A^2*B^2 - 2*A*B*A*B - B^2*A^2 + 2*B*A*B*A)
  assert p1 == A*x + B*y + cxy*x*y + cxy2*x*y^2 + cx2y*x^2*y + cx2y2*x^2*y^2

  

if __name__ == '__main__':
  for f in globals().keys():
    if f.startswith("test_"):
      print '%-20s' % f,
      t0 = time()
      globals()[f]()
      t1 = time()
      print '%.2f' %(t1-t0)

